import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from xgboost import XGBRegressor
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
import matplotlib
matplotlib.use("TkAgg")
warnings.filterwarnings('ignore')

# Configuration Constants
FILE_PATH_POWER = r'D:\study\数学建模\暑假培训\练习题1\第一道练习题\新数据\电站1\电站1发电数据_处理后.xlsx'
FILE_PATH_WEATHER = r'D:\study\数学建模\暑假培训\练习题1\第一道练习题\新数据\电站1\电站 1 天气数据处理后.xlsx'
FILE_PATH_IRRADIANCE = r'D:\study\数学建模\暑假培训\练习题1\第一道练习题\新数据\电站1\电站1环境检测仪数据_填充后_截尾处理后.xlsx'

RESAMPLE_FREQ = '15min'
CLEANING_DATE_MODEL_TRAIN = pd.to_datetime('2024-12-05')
NIGHT_START_HOUR = 22
NIGHT_END_HOUR = 6

# Model Parameters
TEMP_COEFFICIENT = 0.004
REF_TEMP = 25
LAG_HOURS = 1
ROLLING_HOURS = 2

# Alert Thresholds
CLEANING_DROP_THRESHOLD_PCT = -15
CLEANING_DROP_THRESHOLD_ABS = -8
THRESHOLD_MILD_WARNING = 5
THRESHOLD_URGENT_WARNING = 8
CONSECUTIVE_DAYS_MILD = 3
CONSECUTIVE_DAYS_URGENT = 1
TREND_SLOPE_THRESHOLD = 0.5
LOSS_PCT_THRESHOLD = 3

# Cost Analysis Parameters
INSTALLED_CAPACITY_KW = 4998.30
CLEANING_COST_PER_KW = 2
AVERAGE_ELECTRICITY_PRICE_YUAN_PER_KWH = 2
ROI_THRESHOLD = 1.5

def load_data():
    """Load data from Excel files with error handling"""
    try:
        df_power = pd.read_excel(FILE_PATH_POWER)
        df_weather = pd.read_excel(FILE_PATH_WEATHER)
        df_irradiance = pd.read_excel(FILE_PATH_IRRADIANCE)
        print("数据加载成功！(从Excel文件)")
        return df_power, df_weather, df_irradiance
    except FileNotFoundError as e:
        print(f"错误：文件未找到。请检查文件路径和文件名是否正确。\n错误信息: {e}")
        print("将使用模拟数据进行演示。")
        return None, None, None

def prepare_datetime_index(df, time_col='时间'):
    """Convert time column to datetime and set as index"""
    df[time_col] = pd.to_datetime(df[time_col])
    return df.set_index(time_col).sort_index()

def create_time_series_index(df, freq='5min'):
    """Create complete time series index with specified frequency"""
    if df.empty:
        raise ValueError("df数据为空，无法设置频率或进行后续处理。")
    
    min_ts = df.index.min()
    max_ts = df.index.max()
    full_time_idx = pd.date_range(start=min_ts.floor('D'), end=max_ts.ceil('D'), freq=freq)
    return df.reindex(full_time_idx, fill_value=0)

def calculate_instantaneous_power(df_power):
    """Calculate instantaneous power generation with night-time logic"""
    df_power['瞬时发电量kwh'] = np.nan
    valid_days = df_power.index.normalize().unique()
    
    for day in valid_days:
        day_start = pd.Timestamp(day).normalize()
        day_end = day_start + pd.Timedelta(days=1) - pd.Timedelta(minutes=5)
        
        if not df_power.index.to_series().between(day_start, day_end, inclusive='both').any():
            continue
            
        day_df = df_power.loc[day_start:day_end].copy()
        if day_df.empty:
            continue
            
        day_start_time = pd.to_datetime(f'{NIGHT_END_HOUR:02}:00').time()
        day_end_time = pd.to_datetime(f'{NIGHT_START_HOUR-1:02}:59').time()
        
        day_df_daytime = day_df.between_time(day_start_time, day_end_time)
        
        if not day_df_daytime.empty:
            day_df_daytime['瞬时发电量kwh_temp'] = day_df_daytime['当日累计发电量kwh'].diff()
            
            first_valid_idx = day_df_daytime['当日累计发电量kwh'].gt(0).idxmax() if day_df_daytime['当日累计发电量kwh'].gt(0).any() else None
            if first_valid_idx is not None and first_valid_idx in day_df_daytime.index:
                day_df_daytime.loc[first_valid_idx, '瞬时发电量kwh_temp'] = day_df_daytime.loc[first_valid_idx, '当日累计发电量kwh']
            else:
                day_df_daytime['瞬时发电量kwh_temp'] = 0
                
            df_power.loc[day_df_daytime.index, '瞬时发电量kwh'] = day_df_daytime['瞬时发电量kwh_temp']
        
        # Set night-time values to 0
        night_part1_start = day_start + pd.Timedelta(hours=NIGHT_START_HOUR)
        night_part1_end = day_start + pd.Timedelta(hours=23, minutes=59, seconds=59)
        night_part2_start = day_start + pd.Timedelta(days=1)
        night_part2_end = day_start + pd.Timedelta(days=1, hours=NIGHT_END_HOUR-1, minutes=59, seconds=59)
        
        night_mask = (df_power.index.to_series().between(night_part1_start, night_part1_end, inclusive='both') |
                     df_power.index.to_series().between(night_part2_start, night_part2_end, inclusive='both'))
        df_power.loc[night_mask, '瞬时发电量kwh'] = 0
    
    df_power['瞬时发电量kwh'] = np.maximum(df_power['瞬时发电量kwh'], 0).fillna(0)
    return df_power

def resample_data(df_power, df_irradiance, df_weather):
    """Resample all dataframes to consistent frequency"""
    df_power_resampled = df_power['瞬时发电量kwh'].resample(RESAMPLE_FREQ).sum().to_frame()
    df_power_resampled.columns = ['实际发电量kwh']
    
    df_irradiance_resampled = df_irradiance['辐照强度w/m2'].resample(RESAMPLE_FREQ).mean().to_frame()
    
    df_weather_numeric = df_weather[['当前温度', '最高温度', '最低温度', '风速', '湿度']].resample(RESAMPLE_FREQ).mean()
    df_weather_categorical = df_weather[['天气', '风向', '日出时间', '日落时间']].resample(RESAMPLE_FREQ).ffill()
    df_weather_resampled = pd.concat([df_weather_numeric, df_weather_categorical], axis=1)
    
    return df_power_resampled, df_irradiance_resampled, df_weather_resampled

def merge_and_fill_data(df_power, df_irradiance, df_weather):
    """Merge all data and fill missing values"""
    df_merged = pd.concat([df_power, df_irradiance, df_weather], axis=1)
    
    # Set frequency if missing
    if df_merged.index.freq is None:
        inferred_freq = pd.infer_freq(df_merged.index)
        if inferred_freq is None:
            try:
                df_merged.index.freq = RESAMPLE_FREQ
            except ValueError:
                pass
        else:
            df_merged.index.freq = inferred_freq
    
    # Fill numeric columns with interpolation
    numeric_cols = ['辐照强度w/m2', '当前温度', '最高温度', '最低温度', '风速', '湿度']
    for col in numeric_cols:
        if col in df_merged.columns:
            df_merged[col] = df_merged[col].interpolate(method='linear', limit_direction='both')
    
    # Fill categorical columns
    categorical_cols = ['天气', '风向', '日出时间', '日落时间']
    for col in categorical_cols:
        if col in df_merged.columns:
            df_merged[col] = df_merged[col].fillna(method='ffill').fillna(method='bfill')
    
    df_merged['实际发电量kwh'] = df_merged['实际发电量kwh'].fillna(0)
    df_merged['日出时间'] = pd.to_datetime(df_merged['日出时间'])
    df_merged['日落时间'] = pd.to_datetime(df_merged['日落时间'])
    
    return df_merged

def add_time_features(df):
    """Add time-based features"""
    df['小时'] = df.index.hour
    df['一年中的第几天'] = df.index.dayofyear
    df['星期几'] = df.index.dayofweek
    df['月份'] = df.index.month
    return df

def add_interaction_features(df):
    """Add interaction and lag features"""
    df['辐照_当前温度_交互'] = df['辐照强度w/m2'] * df['当前温度']
    df['辐照_效率修正'] = df['辐照强度w/m2'] * (1 - TEMP_COEFFICIENT * (df['当前温度'] - REF_TEMP))
    df['辐照_效率修正'] = np.maximum(df['辐照_效率修正'], 0)
    
    # Lag features
    lag_intervals = int(pd.Timedelta(f'{LAG_HOURS}H') / pd.Timedelta(RESAMPLE_FREQ))
    df['辐照_lag_1h'] = df['辐照强度w/m2'].shift(lag_intervals)
    df['当前温度_lag_1h'] = df['当前温度'].shift(lag_intervals)
    df['风速_lag_1h'] = df['风速'].shift(lag_intervals)
    
    # Rolling features
    rolling_intervals = int(pd.Timedelta(f'{ROLLING_HOURS}H') / pd.Timedelta(RESAMPLE_FREQ))
    df['辐照_rolling_mean_2h'] = df['辐照强度w/m2'].rolling(window=rolling_intervals, min_periods=1).mean()
    
    return df

def infer_cleaning_dates(df_daily_avg_dust_impact):
    """Infer historical cleaning dates from dust impact data"""
    if df_daily_avg_dust_impact.empty:
        return []
    
    daily_change = df_daily_avg_dust_impact.diff().fillna(0)
    cleaning_dates = []
    
    for i in range(1, len(df_daily_avg_dust_impact)):
        current_date = df_daily_avg_dust_impact.index[i]
        prev_impact = df_daily_avg_dust_impact.iloc[i-1]
        current_impact = df_daily_avg_dust_impact.iloc[i]
        
        if (prev_impact > 3 and 
            (daily_change.iloc[i] < CLEANING_DROP_THRESHOLD_ABS or 
             (prev_impact != 0 and (current_impact - prev_impact) / prev_impact * 100 < CLEANING_DROP_THRESHOLD_PCT))):
            cleaning_dates.append(current_date.normalize())
    
    return sorted(list(set(cleaning_dates)))

def get_cleaning_alert(date, df_daily_avg_dust_impact, df_full_data):
    """Generate cleaning alerts for a specific date"""
    alerts = []
    current_date = pd.to_datetime(date).normalize()
    
    if current_date not in df_daily_avg_dust_impact.index:
        alerts.append(f"警告：{current_date.strftime('%Y-%m-%d')} 无日均积灰指数数据。")
        return alerts
    
    current_impact = df_daily_avg_dust_impact.loc[current_date]
    alerts.append(f"在 {current_date.strftime('%Y-%m-%d')}，日均积灰影响指数为 {current_impact:.2f}%。")
    
    # Urgent warning check
    start_date_urgent = current_date - pd.Timedelta(days=CONSECUTIVE_DAYS_URGENT - 1)
    recent_impact_urgent = df_daily_avg_dust_impact.loc[start_date_urgent:current_date]
    
    if (not recent_impact_urgent.empty and len(recent_impact_urgent) == CONSECUTIVE_DAYS_URGENT and
        (recent_impact_urgent > THRESHOLD_URGENT_WARNING).all()):
        alerts.append(f"！！！紧急预警！！！：连续 {CONSECUTIVE_DAYS_URGENT} 天日均积灰影响指数超过 {THRESHOLD_URGENT_WARNING}%。建议立即安排清洗！")
    
    # Trend analysis
    trend_start = current_date - pd.Timedelta(days=7)
    trend_data = df_daily_avg_dust_impact.loc[trend_start:current_date].dropna()
    
    if len(trend_data) >= 3 and len(trend_data.unique()) > 1:
        x = np.arange(len(trend_data))
        y = trend_data.values
        slope, _ = np.polyfit(x, y, 1)
        if slope > TREND_SLOPE_THRESHOLD:
            alerts.append(f"趋势预警：积灰影响指数呈上升趋势 (日均增长约 {slope:.2f}%)。")
    
    # Mild warning check
    start_date_mild = current_date - pd.Timedelta(days=CONSECUTIVE_DAYS_MILD - 1)
    recent_impact_mild = df_daily_avg_dust_impact.loc[start_date_mild:current_date]
    
    if (not recent_impact_mild.empty and len(recent_impact_mild) == CONSECUTIVE_DAYS_MILD and
        (recent_impact_mild > THRESHOLD_MILD_WARNING).all()):
        alerts.append(f"清洗建议：连续 {CONSECUTIVE_DAYS_MILD} 天日均积灰影响指数超过 {THRESHOLD_MILD_WARNING}%。建议考虑安排清洗。")
    elif current_impact <= THRESHOLD_MILD_WARNING:
        alerts.append(f"目前日均积灰影响指数低于 {THRESHOLD_MILD_WARNING}%，情况良好。")
    
    return alerts

def get_cost_benefit_analysis(date, df_daily_avg_dust_impact, df_full_data):
    """Perform cost-benefit analysis for cleaning decision"""
    alerts = []
    current_date = pd.to_datetime(date).normalize()
    
    total_cleaning_cost = INSTALLED_CAPACITY_KW * CLEANING_COST_PER_KW
    alerts.append(f"总清洗成本: {total_cleaning_cost:.2f} 元")
    
    if current_date not in df_daily_avg_dust_impact.index:
        alerts.append(f"警告：{current_date.strftime('%Y-%m-%d')} 无数据。")
        return alerts
    
    current_impact = df_daily_avg_dust_impact.loc[current_date]
    if current_impact < 6:
        alerts.append(f"积灰影响指数 ({current_impact:.2f}%) 较低，暂不建议清洗。")
        return alerts
    
    # Calculate potential benefit
    projected_days = 7
    avg_daily_power = 100  # Simplified assumption
    projected_loss_revenue = 0
    
    for i in range(1, projected_days + 1):
        projected_impact = min(current_impact + i * 1, 100)
        daily_loss = avg_daily_power * (projected_impact / 100)
        projected_loss_revenue += daily_loss * AVERAGE_ELECTRICITY_PRICE_YUAN_PER_KWH
    
    alerts.append(f"未来 {projected_days} 天预估损失: {projected_loss_revenue:.2f} 元")
    
    if total_cleaning_cost > 0:
        roi = projected_loss_revenue / total_cleaning_cost
        alerts.append(f"投资回报率: {roi:.2f}")
        
        if roi >= ROI_THRESHOLD:
            alerts.append("***建议清洗***：效益大于成本")
        else:
            alerts.append("***暂不建议清洗***：效益低于成本")
    
    return alerts

def visualize_results(df, df_daily_avg_dust_impact, historical_cleaning_dates):
    """Create visualization plots"""
    plt.style.use('seaborn-v0_8-darkgrid')
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    
    # Plot 1: Power generation comparison
    plt.figure(figsize=(20, 10))
    plt.plot(df.index, df['实际发电量kwh'], label='实际发电量 (kWh)', alpha=0.8)
    plt.plot(df.index, df['预测发电量kwh'], label='预测理论发电量 (kWh)', linestyle='--', alpha=0.8)
    plt.fill_between(df.index, df['实际发电量kwh'], df['预测发电量kwh'],
                     where=(df['预测发电量kwh'] > df['实际发电量kwh']),
                     color='red', alpha=0.2, label='积灰损失 (kWh)')
    
    plt.axvline(x=CLEANING_DATE_MODEL_TRAIN, color='green', linestyle=':', linewidth=2,
                label=f'模型训练基准清洗日期: {CLEANING_DATE_MODEL_TRAIN.strftime("%Y-%m-%d")}')
    
    for i, clean_date in enumerate(historical_cleaning_dates):
        plt.axvline(x=clean_date, color='blue', linestyle='--', linewidth=1, alpha=0.7,
                    label='历史清洗日期' if i == 0 else "")
    
    plt.title('光伏电站发电量趋势与积灰损失')
    plt.xlabel('时间')
    plt.ylabel('发电量 (kWh)')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()
    
    # Plot 2: Dust impact index
    plt.figure(figsize=(20, 8))
    plt.plot(df_daily_avg_dust_impact.index, df_daily_avg_dust_impact.values,
             label='日均积灰影响指数 (%)', color='orange', alpha=0.8)
    
    for i, clean_date in enumerate(historical_cleaning_dates):
        plt.axvline(x=clean_date, color='blue', linestyle='--', linewidth=1, alpha=0.7,
                    label='历史清洗日期' if i == 0 else "")
    
    plt.axvline(x=CLEANING_DATE_MODEL_TRAIN, color='green', linestyle=':', linewidth=2,
                label=f'模型训练基准清洗日期: {CLEANING_DATE_MODEL_TRAIN.strftime("%Y-%m-%d")}')
    
    plt.title('积灰影响指数随时间变化 (日均)')
    plt.xlabel('日期')
    plt.ylabel('积灰影响指数 (%)')
    plt.axhline(y=5, color='gray', linestyle='--', label='清洗建议阈值 (5%)')
    plt.axhline(y=8, color='red', linestyle='--', label='紧急清洗阈值 (8%)')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

def main():
    """Main execution function"""
    # Load and prepare data
    df_power, df_weather, df_irradiance = load_data()
    if df_power is None:
        return
    
    # Prepare time indices
    df_power = prepare_datetime_index(df_power)
    df_weather = prepare_datetime_index(df_weather)
    df_irradiance = prepare_datetime_index(df_irradiance)
    
    # Create complete time series
    df_power = create_time_series_index(df_power)
    df_power['当日累计发电量kwh'] = df_power['当日累计发电量kwh'].replace(0, np.nan).fillna(method='ffill').fillna(method='bfill').fillna(0)
    
    # Calculate instantaneous power
    df_power = calculate_instantaneous_power(df_power)
    
    # Resample data
    df_power_resampled, df_irradiance_resampled, df_weather_resampled = resample_data(df_power, df_irradiance, df_weather)
    
    # Merge and process data
    df_merged = merge_and_fill_data(df_power_resampled, df_irradiance_resampled, df_weather_resampled)
    df = df_merged.copy()
    
    # Feature engineering
    df = add_time_features(df)
    df['is_daytime'] = (df.index.time >= df['日出时间'].dt.time) & (df.index.time <= df['日落时间'].dt.time)
    
    # Set night-time irradiance to 0
    is_physical_night = (df.index.hour >= df['日落时间'].dt.hour) | (df.index.hour < df['日出时间'].dt.hour)
    df.loc[is_physical_night, '辐照强度w/m2'] = 0
    
    df_daytime = df[df['is_daytime']].copy()
    if df_daytime.empty:
        print("警告：没有有效的白天数据")
        return
    
    # Set frequency for daytime data
    if df_daytime.index.freq is None:
        min_ts = df_daytime.index.min()
        max_ts = df_daytime.index.max()
        full_time_idx = pd.date_range(start=min_ts, end=max_ts, freq=RESAMPLE_FREQ)
        df_daytime = df_daytime.reindex(full_time_idx)
        df_daytime.index.freq = RESAMPLE_FREQ
    
    df_daytime = add_interaction_features(df_daytime)
    
    # Categorical encoding
    df_daytime = pd.get_dummies(df_daytime, columns=['天气', '风向'], drop_first=True)
    
    # Add wind direction interactions
    wind_direction_cols = [col for col in df_daytime.columns if col.startswith('风向_')]
    for col in wind_direction_cols:
        wind_direction_name = col.replace('风向_', '')
        df_daytime[f'风速_{wind_direction_name}_交互'] = df_daytime['风速'] * df_daytime[col]
        df_daytime[f'湿度_{wind_direction_name}_交互'] = df_daytime['湿度'] * df_daytime[col]
    
    df_daytime.dropna(inplace=True)
    if df_daytime.empty:
        print("警告：特征工程后无有效数据")
        return
    
    # Train model
    clean_start = CLEANING_DATE_MODEL_TRAIN + pd.Timedelta(days=7)
    clean_end = CLEANING_DATE_MODEL_TRAIN + pd.Timedelta(days=14)
    df_clean = df_daytime[(df_daytime.index >= clean_start) & (df_daytime.index <= clean_end)].copy()
    
    features = [
        '辐照强度w/m2', '当前温度', '最高温度', '最低温度', '风速', '湿度',
        '小时', '一年中的第几天', '星期几', '月份',
        '辐照_当前温度_交互', '辐照_效率修正',
        '辐照_lag_1h', '当前温度_lag_1h', '风速_lag_1h', '辐照_rolling_mean_2h'
    ]
    features.extend([col for col in df_daytime.columns if col.startswith(('天气_', '风向_', '风速_', '湿度_')) and col.endswith('_交互')])
    
    common_features = list(set(features) & set(df_daytime.columns))
    X_train = df_clean[common_features]
    y_train = df_clean['实际发电量kwh']
    
    if X_train.empty:
        raise ValueError("没有足够的干净数据用于训练模型")
    
    # Train XGBoost model
    xgb_model = XGBRegressor(
        random_state=42, n_estimators=300, learning_rate=0.05,
        max_depth=7, subsample=0.8, colsample_bytree=0.8, n_jobs=-1
    )
    xgb_model.fit(X_train, y_train)
    
    # Make predictions
    X_pred = df_daytime[common_features]
    df_daytime['预测发电量kwh'] = np.maximum(xgb_model.predict(X_pred), 0)
    
    df['预测发电量kwh'] = 0.0
    df.loc[df_daytime.index, '预测发电量kwh'] = df_daytime['预测发电量kwh']
    
    # Calculate dust impact
    df['发电量损失kwh'] = np.maximum(df['预测发电量kwh'] - df['实际发电量kwh'], 0)
    df['积灰影响指数'] = (df['发电量损失kwh'] / df['预测发电量kwh']).fillna(0) * 100
    df.loc[df['预测发电量kwh'] < 0.1, '积灰影响指数'] = 0
    
    df_daily_avg_dust_impact = df[df['is_daytime'] & (df['预测发电量kwh'] > 0.1)]['积灰影响指数'].resample('D').mean().dropna()
    
    # Infer cleaning dates
    historical_cleaning_dates = infer_cleaning_dates(df_daily_avg_dust_impact)
    
    # Results summary
    print("\n--- 关键结果摘要 ---")
    print(f"数据总时间范围: {df.index.min()} 到 {df.index.max()}")
    print(f"总实际发电量: {df['实际发电量kwh'].sum():.2f} kWh")
    print(f"总预测理论发电量: {df['预测发电量kwh'].sum():.2f} kWh")
    print(f"总积灰造成的发电量损失: {df['发电量损失kwh'].sum():.2f} kWh")
    print(f"平均积灰影响指数: {df_daily_avg_dust_impact.mean():.2f}%")
    
    if historical_cleaning_dates:
        print(f"推断出 {len(historical_cleaning_dates)} 条历史清洗事件")
    
    # Visualizations
    visualize_results(df, df_daily_avg_dust_impact, historical_cleaning_dates)
    
    # Alert analysis for recent dates
    print("\n--- 实时清洗预警规则 ---")
    evaluable_dates = df_daily_avg_dust_impact.index.normalize().unique().sort_values()
    recent_dates = evaluable_dates[-15:] if len(evaluable_dates) >= 15 else evaluable_dates
    
    for current_date in recent_dates:
        if current_date in df_daily_avg_dust_impact.index:
            print(f"\n--- {current_date.strftime('%Y-%m-%d')} 清洗预警评估 ---")
            alerts = get_cleaning_alert(current_date, df_daily_avg_dust_impact, df)
            for alert in alerts:
                print(f"- {alert}")
    
    # Cost-benefit analysis
    print("\n--- 基于成本效益的清洗决策 ---")
    for current_date in recent_dates:
        if current_date in df_daily_avg_dust_impact.index:
            print(f"\n--- {current_date.strftime('%Y-%m-%d')} 成本效益评估 ---")
            cost_alerts = get_cost_benefit_analysis(current_date, df_daily_avg_dust_impact, df)
            for alert in cost_alerts:
                print(f"- {alert}")

if __name__ == "__main__":
    main()